void rattle2(celldata_t *cell, int i){
  int ig0 = i;
  int ig1 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[0]];
  real rmass0 = cell->rmass[ig0];
  real rmass1 = cell->rmass[ig1];
  vec<real> r01, vp01;
  vecsubv(r01, cell->x[ig1], cell->x[ig0]);
  vecsubv(vp01, cell->shake_vp[ig1], cell->shake_vp[ig0]);
  real a00 = vecdot(r01, r01)*(+rmass0+rmass1);
  real invmatdet = 1/(a00);
  real a00inv = (1) * invmatdet;
  real c0 = -vecdot(vp01, r01);
  real l0 = a00inv*c0;
  vecscaleaddv(cell->v[ig0], cell->v[ig0], r01, 1, -l0 * rmass0);
  vecscaleaddv(cell->v[ig1], cell->v[ig1], r01, 1, +l0 * rmass1);
}
void rattle3(celldata_t *cell, int i){
  int ig0 = i;
  int ig1 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[0]];
  int ig2 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[1]];
  real rmass0 = cell->rmass[ig0];
  real rmass1 = cell->rmass[ig1];
  real rmass2 = cell->rmass[ig2];
  vec<real> r01, vp01, r02, vp02;
  vecsubv(r01, cell->x[ig1], cell->x[ig0]);
  vecsubv(r02, cell->x[ig2], cell->x[ig0]);
  vecsubv(vp01, cell->shake_vp[ig1], cell->shake_vp[ig0]);
  vecsubv(vp02, cell->shake_vp[ig2], cell->shake_vp[ig0]);
  real a00 = vecdot(r01, r01)*(+rmass0+rmass1);
  real a01 = vecdot(r01, r02)*(+rmass0);
  real a10 = a01;
  real a11 = vecdot(r02, r02)*(+rmass0+rmass2);
  real invmatdet = 1/(a00*a11 - (a01)*(a01));
  real a00inv = (a11) * invmatdet;
  real a01inv = (-a01) * invmatdet;
  real a11inv = (a00) * invmatdet;
  real c0 = -vecdot(vp01, r01);
  real c1 = -vecdot(vp02, r02);
  real l0 = a00inv*c0 + a01inv*c1;
  real l1 = a01inv*c0 + a11inv*c1;
  vecscaleaddv(cell->v[ig0], cell->v[ig0], r01, 1, -l0 * rmass0);
  vecscaleaddv(cell->v[ig0], cell->v[ig0], r02, 1, -l1 * rmass0);
  vecscaleaddv(cell->v[ig1], cell->v[ig1], r01, 1, +l0 * rmass1);
  vecscaleaddv(cell->v[ig2], cell->v[ig2], r02, 1, +l1 * rmass2);
}
void rattle3angle(celldata_t *cell, int i){
  int ig0 = i;
  int ig1 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[0]];
  int ig2 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[1]];
  real rmass0 = cell->rmass[ig0];
  real rmass1 = cell->rmass[ig1];
  real rmass2 = cell->rmass[ig2];
  vec<real> r01, vp01, r02, vp02, r12, vp12;
  vecsubv(r01, cell->x[ig1], cell->x[ig0]);
  vecsubv(r02, cell->x[ig2], cell->x[ig0]);
  vecsubv(r12, cell->x[ig2], cell->x[ig1]);
  vecsubv(vp01, cell->shake_vp[ig1], cell->shake_vp[ig0]);
  vecsubv(vp02, cell->shake_vp[ig2], cell->shake_vp[ig0]);
  vecsubv(vp12, cell->shake_vp[ig2], cell->shake_vp[ig1]);
  real a00 = vecdot(r01, r01)*(+rmass0+rmass1);
  real a01 = vecdot(r01, r02)*(+rmass0);
  real a02 = vecdot(r01, r12)*(-rmass1);
  real a10 = a01;
  real a11 = vecdot(r02, r02)*(+rmass0+rmass2);
  real a12 = vecdot(r02, r12)*(+rmass2);
  real a20 = a02;
  real a21 = a12;
  real a22 = vecdot(r12, r12)*(+rmass1+rmass2);
  real invmatdet = 1/(a00*a11*a22 - a00*(a12)*(a12) - (a01)*(a01)*a22 + 2*a01*a02*a12 - (a02)*(a02)*a11);
  real a00inv = (a11*a22 - (a12)*(a12)) * invmatdet;
  real a01inv = (-a01*a22 + a02*a12) * invmatdet;
  real a02inv = (a01*a12 - a02*a11) * invmatdet;
  real a11inv = (a00*a22 - (a02)*(a02)) * invmatdet;
  real a12inv = (-a00*a12 + a01*a02) * invmatdet;
  real a22inv = (a00*a11 - (a01)*(a01)) * invmatdet;
  real c0 = -vecdot(vp01, r01);
  real c1 = -vecdot(vp02, r02);
  real c2 = -vecdot(vp12, r12);
  real l0 = a00inv*c0 + a01inv*c1 + a02inv*c2;
  real l1 = a01inv*c0 + a11inv*c1 + a12inv*c2;
  real l2 = a02inv*c0 + a12inv*c1 + a22inv*c2;
  vecscaleaddv(cell->v[ig0], cell->v[ig0], r01, 1, -l0 * rmass0);
  vecscaleaddv(cell->v[ig0], cell->v[ig0], r02, 1, -l1 * rmass0);
  vecscaleaddv(cell->v[ig1], cell->v[ig1], r01, 1, +l0 * rmass1);
  vecscaleaddv(cell->v[ig1], cell->v[ig1], r12, 1, -l2 * rmass1);
  vecscaleaddv(cell->v[ig2], cell->v[ig2], r02, 1, +l1 * rmass2);
  vecscaleaddv(cell->v[ig2], cell->v[ig2], r12, 1, +l2 * rmass2);
}
void rattle4(celldata_t *cell, int i){
  int ig0 = i;
  int ig1 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[0]];
  int ig2 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[1]];
  int ig3 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[2]];
  real rmass0 = cell->rmass[ig0];
  real rmass1 = cell->rmass[ig1];
  real rmass2 = cell->rmass[ig2];
  real rmass3 = cell->rmass[ig3];
  vec<real> r01, vp01, r02, vp02, r03, vp03;
  vecsubv(r01, cell->x[ig1], cell->x[ig0]);
  vecsubv(r02, cell->x[ig2], cell->x[ig0]);
  vecsubv(r03, cell->x[ig3], cell->x[ig0]);
  vecsubv(vp01, cell->shake_vp[ig1], cell->shake_vp[ig0]);
  vecsubv(vp02, cell->shake_vp[ig2], cell->shake_vp[ig0]);
  vecsubv(vp03, cell->shake_vp[ig3], cell->shake_vp[ig0]);
  real a00 = vecdot(r01, r01)*(+rmass0+rmass1);
  real a01 = vecdot(r01, r02)*(+rmass0);
  real a02 = vecdot(r01, r03)*(+rmass0);
  real a10 = a01;
  real a11 = vecdot(r02, r02)*(+rmass0+rmass2);
  real a12 = vecdot(r02, r03)*(+rmass0);
  real a20 = a02;
  real a21 = a12;
  real a22 = vecdot(r03, r03)*(+rmass0+rmass3);
  real invmatdet = 1/(a00*a11*a22 - a00*(a12)*(a12) - (a01)*(a01)*a22 + 2*a01*a02*a12 - (a02)*(a02)*a11);
  real a00inv = (a11*a22 - (a12)*(a12)) * invmatdet;
  real a01inv = (-a01*a22 + a02*a12) * invmatdet;
  real a02inv = (a01*a12 - a02*a11) * invmatdet;
  real a11inv = (a00*a22 - (a02)*(a02)) * invmatdet;
  real a12inv = (-a00*a12 + a01*a02) * invmatdet;
  real a22inv = (a00*a11 - (a01)*(a01)) * invmatdet;
  real c0 = -vecdot(vp01, r01);
  real c1 = -vecdot(vp02, r02);
  real c2 = -vecdot(vp03, r03);
  real l0 = a00inv*c0 + a01inv*c1 + a02inv*c2;
  real l1 = a01inv*c0 + a11inv*c1 + a12inv*c2;
  real l2 = a02inv*c0 + a12inv*c1 + a22inv*c2;
  vecscaleaddv(cell->v[ig0], cell->v[ig0], r01, 1, -l0 * rmass0);
  vecscaleaddv(cell->v[ig0], cell->v[ig0], r02, 1, -l1 * rmass0);
  vecscaleaddv(cell->v[ig0], cell->v[ig0], r03, 1, -l2 * rmass0);
  vecscaleaddv(cell->v[ig1], cell->v[ig1], r01, 1, +l0 * rmass1);
  vecscaleaddv(cell->v[ig2], cell->v[ig2], r02, 1, +l1 * rmass2);
  vecscaleaddv(cell->v[ig3], cell->v[ig3], r03, 1, +l2 * rmass3);
}
